Chebyshev approach to quantum systems coupled to a bath 
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We propose a new concept for the dynamics of a quantum bath, the Chebyshev space, and a new 
method based on this concept, the Chebyshev space method. The Chebyshev space is an abstract 
vector space that exactly represents the fermionic or bosonic bath degrees of freedom, without a 
discretization of the bath density of states. Relying on Chebyshev expansions the Chebyshev space 
representation of a bath has very favorable properties with respect to extremely precise and efficient 
calculations of groundstate properties, static and dynamical correlations, and time-evolution for a 
great variety of quantum systems. The aim of the present work is to introduce the Chebyshev 
space in detail and to demonstrate the capabilities of the Chebyshev space method. Although 
the central idea is derived in full generality the focus is on model systems coupled to fermionic 
baths. In particular we address quantum impurity problems, such as an impurity in a host or a 
bosonic impurity with a static barrier, and the motion of a wave packet on a chain coupled to leads. 
For the bosonic impurity, the phase transition from a delocalized electron to a localized polaron 
in arbitrary dimension is detected. For the wave packet on a chain, we show how the Chebyshev 
space method implements different boundary conditions, including transparent boundary conditions 
replacing infinite leads. Furthermore the self-consistent solution of the Holstein model in infinite 
dimension is calculated. With the examples we demonstrate how highly accurate results for system 
energies, correlation and spectral functions, and time-dependence of observables are obtained with 
modest computational effort. 

PACS numbers: 71.27.+a,73.21.-b,73.23.-b 



I. INTRODUCTION 



The calculation of spectral or dynamical properties of 
quantum systems, expressed through spectral functions 
or captured in the time evolution of a wave function, 
is one of the most important and most promising ap- 
plications of modern numerical techniques in theoretical 
physics or chemistry. For many purposes, approxima- 
tion free techniques that allow to calculate numerically 
exact results for arbitrary Hamiltonians on large Hilbcrt 
spaces are of interest. One of the most powerful tools in 
this context are techniques based on Chebyshev expan- 
sions, like the kernel polynomial method (KPM)i*2i2iii^ ) 
which yield results of high accuracy with modest compu- 
tational effort. Chebyshev techniques often outperform 
the Lanczos algorithm with respect to accuracy and ef- 
ficiency, and are applicable to problems that are beyond 
the reach of matrix diagonalization. In this article we 
propose an extension of Chebyshev techniques that con- 
siderably enlarges their field of applications. Possible 
new applications we have in mind include (i) damping 
and decoherence in quantum systems coupled to an envi- 
ronment like a bosonic heat bath, (ii) transport through 
quantum systems coupled to fermionic baths, (iii) the so- 
lution of quantum impurity models, (iv) non-equilibrium 
dynamics of mesoscopic devices coupled to leads, and 
in general, (v) the calculation of static and dynamical 
correlation functions or time propagation in these phys- 
ical situations, (vi) the combination with diagrammatic 
Green function techniques, (vii) the treatment of degrees 
of freedom with non-trivial dynamics like phonons with 
dispersion. 

We concentrate here on a situation for which the de- 



velopment of ideas is particularly clear: A quantum sys- 
tem coupled to a fermionic bath. The bath serves as a 
reservoir for fermions, which can hop from the bath to 
the quantum system and back. The standard example is 
that of a mesoscopic system contacted with leads For a 
single quantum dot, the appropriate model is the famous 
Anderson model for the Kondo effect^, which describes 
an impurity site with Coulomb interaction embedded in a 
host of non-interacting electrons. This model also arises 
in dynamical mean-field theory (DMFT) where the so- 
lution of an Anderson model in dependence on a host 
spectral function is a central issued 

In all these cases the quantum system allows for many- 
particle interactions, which give rise to non-trivial corre- 
lations. For their description one has to rely on several 
correlation functions whose calculation requires the full 
apparatus of many-body theory. In the numerical cal- 
culation we have to deal with complicated objects like 
a Fock space of many fermions. Due to the coupling 
between the quantum system and the bath, further cor- 
relations between the quantum system and the bath de- 
velop. As a consequence one cannot separate the degrees 
of freedom of the quantum system and of the bath, but 
must treat their evolution in parallel. However, the bath 
consists of non-interacting fermions, and is entirely de- 
scribed by single-particle spectral functions, as no cor- 
relations exist in absence of coupling to the quantum 
system. Plainly spoken, we do not need to describe all 
details of the bath, but only how its presence influences 
the dynamics of the quantum system. It is important 
to exploit this simplification to make a calculation pos- 
sible. Trivial as this point may seem, it is not easily 
incorporated into a numerical calculation. The numeri- 
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FIG. 1: Different system-bath geometries used for examples in 
this article. On the left, a single site is embedded in a host (cf. 
Eq. (HE}), e.g. an impurity in a lattice. In the middle, a chain 
is coupled to a bath (cf. Eq. ([28}), allowing an electron to 
move along the chain and to hop to the bath and return. On 
the right, a single site with local bosonic degrees of freedom 
is coupled to a bath (cf. Eq. ((35}). An electron can excite 
bosons at the site, and hop to the bath while the bosons are 
left behind. The last two cases fit into the scheme of Eq. ((6} , 
where the coupling between system and bath is given by a 
Hamiltonian Hl- 



cally exact techniques to which Chebyshev techniques be- 
long require an explicit representation of the bath degrees 
of freedom in terms of a Hamiltonian and an associated 
Hilbert space. In the worst case such a representation 
again involves a many-body Fock space, describing the 
bath with the same complexity as the interacting quan- 
tum system. 

The new idea we promote here provides a different so- 
lution which is based on Chebyshev expansions of spec- 
tral functions of the bath. The Chebyshev expansion sup- 
plies the information necessary for an exact calculation of 
correlations function for the quantum system coupled to 
the bath, but avoids the introduction of redundant pa- 
rameters that unnecessarily complicate the calculation. 
Thereby, the simplification of a non-interacting bath is 
exploited while the full coupled dynamics of the quan- 
tum system and the bath is retained. The realization of 
this idea amounts to the construction of an abstract vec- 
tor space, which we call the Chebyshev space (CS). The 
combination of the CS with computational techniques 
based on Chebyshev expansions, in particular KPM, will 
be referred to as the Chebyshev space method (CSM). 

We demonstrate the CSM in this article for a number 
of examples (cf. Fig. [T]), including (i) the calculation of 
spectral functions and groundstate energy for an impu- 
rity in a host in various dimensions, (ii) the description of 
the electron-polaron phase transition for a bosonic impu- 
rity, (iii) the self-consistent solution of the Holstein model 
within DMFT, and (iv) the time propagation of a wave 
packet on a chain with different boundaries conditions. 

We do not include examples for interacting fermions 
at finite density that involve renormalization of the bath 
due to the creation of particle-hole pairs like in the Kondo 
problem. The discussion would exceed a tolerable length 
of this article and obscure the presentation of the new 
ideas. We return to this issue in the conclusion, and 
concentrate here on examples better suited for the intro- 
duction of the technical details of the CS construction. 

The article is organized as follows. In Sec. HH we re- 
mind the reader of Chebyshev expansions for spectral 
functions, and fix some notations used in the remainder. 



In Sec. IIIII we introduce the CS, and explain in Sec. IIVI 
the implementation of the CSM for the example of an 
impurity embedded in a host. In Sec. [V] we derive an 
important property of the CS and discuss its practical 
relevance. Sec. IVII describes how self-consistent calcula- 
tions are performed in CSM, and Sec. IVHl demonstrates 
how the concurrent dynamics of a quantum system and 
a bath, when degrees of freedom evolve in parallel, is 
accounted for. This provides the basis for the study of 
the phase transition for a bosonic impurity in Sec. IVIII1 
and to the self-consistent solution of the Holstein model 
within DMFT in Scc. lVIIIAl In Sec. US we address the 
time evolution of a wave packet on a long chain with dif- 
ferent boundary conditions realized by CSM. We summa- 
rize in Sec.|X] an d point out possible future advancements 
of CSM. In the Appendices, we provide a short account 
of technical aspects of Chebyshev expansions and KPM, 
and the derivation of some mathematical results used in 
the text. 



II. CHEBYSHEV EXPANSIONS AND THE 
KERNEL POLYNOMIAL 

The recurrent theme of this article is the expansion of 
a spectral function in a series of Chebyshev polynomials. 
To make the article self-contained we add an appendix 
on technical issues of Chebyshev expansions and KPM 
(App. [X}. For a more detailed exposition, we refer the 
reader to the recent review Ref. [5|, and references cited 
therein. 

The Chebyshev polynomials T n (x) are defined by the 
two-term recurrence^ 



T (x) = 1 , Ti(x) = x 
T n+1 (x) = 2xT n (x) - T„_i(x) 



(1) 



which is equivalent to T n (x) = cos(n arccosx). The 
Chebyshev polynomials are mutually orthogonal on the 
interval [— 1 , 1] with respect to the scalar product given 
by the weighting function (1 — x 2 ) -1 / 2 , obeying the rela- 
tions 



1 T m (x)T n (x) 
-l ir\/l — x 2 



dx 



1, 

2 u mn i 



o 



(2) 



Eq. ([2]) implies that the Chebyshev polynomials form 
an orthogonal basis for functions defined on the interval 
[—1, 1]. To expand a spectral function, defined as 

A(u>) = {i>\6(u-H)\ip) = — lim Im(^| [c^+i?/-^ 1 \ip) 

(3) 

to some Hamiltonian H and vector \ip), in a series of 
Chebyshev polynomials it must therefore vanish outside 
[—1,1]. This can be achieved by re-scaling the Hamilto- 
nian H as H = pH + q, where p, q is chosen in such a way 
that all eigenvalues of H lie in [—1,1], or equivalently 
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\\H\\ < 1. After a corresponding variable substitution 
x = (lu — q)/p, the Chebyshev expansion for A(u>) reads 



A(w) 



1 



1 



P 1X\f\ 



Mo 



(4) 



for A(uS) on the interval [w m i n , w max ] with cJ m i n = — p + q, 
w max = p + q. Comparing this expansion to the defi- 
nition ([3]), and using the orthogonality relation (J2J), we 
find that the coefficients in this series - the Chebyshev 
moments - are given by 



T n (x)A(px + q)dx = (V>|T n (#)|V>> 



(5) 



We will assume in the following that any Hamiltonian is 
properly scaled to [—1,1]. The introduction of scaling 
factors is straightforward most times. Whenever impor- 
tant, we will explicitly discuss the consequences of scal- 
ing. 

The calculation of the vectors T n (H)\ip) in Eq. can 
be accomplished by means of the two-term recurrence ((T|) . 
Starting with the vector one iteratively obtains the 
(n+l) th vector from the n th and (n— l) th vector. The cal- 
culation of the moments therefore requires (i) to apply H 
to a vector, i.e. to perform matrix-vector multiplication 
if H is explicitly given as a matrix, and (ii) to evaluate a 
scalar product of two vectors. 

In any practical application only a finite number TV 
of moments can be calculated, and A(ui) has to be re- 
constructed from a truncated series, taking the first N 
terms in Eq. ((¥]). The reconstruction of A(u>) from such 
a finite series is the issue of KPM (see App. [A"]) . High 
resolution of A(uS) is already obtained with a fairly small 
number of moments, thus allowing for accurate calcula- 
tions with moderate demands on computational time or 
memory. The resolution of KPM scales linearly with TV. 
It is usually much better than for the Lanczos (Recur- 
sion) Methodii^^, and can be unrestrictedly increased 
without problems 5 . 



III. THE GENERAL SCHEME 

A coupled quantum system and bath is generally de- 
scribed by a Hamiltonian of the forrr&i^ 



H = H S + H L +H l 



(6) 



Hs and Hb denote the Hamiltonian for the quantum sys- 
tem and the bath, respectively, and Hl describes the cou- 
pling between system and bath. Without coupling, for 
Hl = 0, the degrees of freedom of the quantum system 
and the bath evolve independently, i.e. Hs and Hb com- 
mute. We restrict ourselves in this article to situations 
with a single spinless fermion. This allows to introduce 
the CS without a discussion of complications inherent to 
many-fermion physics. The CSM can be extended to fi- 
nite fermion density, or the inclusion of spin degrees of 



freedom, with the same implementation of a CS as given 
in this article. We return to that issue in Sec. |Xj 

Since the bath consists of non-interacting fermions, Hb 
is a bilinear Hamiltonian 



H, 



(7) 



a/3 



where fa are fermionic operators for the bath degrees 
of freedom. The indices a, can denote e.g. an arbitrary 
set of orbitals, sites of a lattice, or, if (Hb)ciP is diagonal, 
the eigenstates of Hb- A change of indices corresponds 
to a unitary transformation of the matrix (-ffs)a/3- The 
system Hamiltonian Hs can be of arbitrary form, involv- 
ing any type of interaction. We do not specify Hs now, 
and discuss examples for various different Hs later. 

More can be said about the coupling term Hl- In 
general, the bath is in contact to one (or a few) site(s) of 
the quantum system, say site with fermionic operator 
Cq . Fermions hop from this site to the bath, and back. 
The appropriate choice for Hl is 

Hl=J2 l Mf a + ftc ) (t a e R) , (8) 

a 

or a sum of terms of this form. It is convenient to intro- 
duce a fermionic operator by 



/2 



(9) 



As X)a a a = 1' the fermionic anticommutator relation 
{d, d^} = 1 is fulfilled. With the introduction of d (t) , H L 
acquires the form of a hopping term between two sites, 



H L = V(cld + Jc ) , V = (£tl) 1 



/2 



(10) 



Depending on the meaning of indices in Eq. l[Tj). d can 
denote a concrete orbital or a site in a lattice, or just an 
abstract linear combination of /-operators. If we think 
of a mesoscopic system contacted to a lead, then Cq^ and 
d't) denote operators for the contact point in the system 
and lead, respectively. 

We define the bath spectral function Ab(uj) as the 
spectral function of the e?-orbital, 



Ab(v) = (va,c\d5(ui — Hb) d'\v&c) 



(11) 



where |vac) is the bath vacuum, i.e. d|vac) = 0. Remem- 
ber that we study situations with a single fermion. For 
finite fermion density, we had to consider both the parti- 
cle and hole part of the spectral function, and to account 
for Pauli blocking. 

After the transformation @, the bath occurs in the 
calculation of correlation functions for the quantum 
system only through the spectral function Ab{w) and 
the coupling strength V. The precise coefficients in 
Eqs. 0, © do not occur. In this sense, all Hamilto- 
nians © with the same Ab{oj), V, but potentially dif- 
ferent Hb, Hl, are equivalent. Often the problem under 
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study is denned in this way, by specifying A B (u) and 
V, without an explicit representation of Hb or Hi, as in 
Eqs. 0, © at hand. 

Given such a problem, how can it be accessed within 
numerically exact techniques like Lanczos or KPM? 
Since these techniques need the Hamiltonian H explicitly 
given, it seems inevitable to use an explicit representa- 
tion of Hb in the form ((7|). One possible way to obtain 
such a representation is to discretize Ab{uj) by a finite 
number of (5-peaks, as 



and 



A B {u) 



]w a 5{uj - e a ) 



(12) 



This approximation translates to Hb — £ a f a fa, and 



d) = J2 a w » 2 fa- A calculation then relies on our abil- 



ity to construct a good approximation (fT2j) . This poses 
certain questions: How to choose the positions e a of 
peaks, how to choose their weights w a , and how many of 
them are needed to approximate Ag(w) sufficiently well? 
When is an approximation to Ab{lo) sufficiently good, 
and what is the precise meaning of in Eq. (TT2"|) ? Is 
there an optimal way to choose e a , w a for a given number 
of peaks? While the result of a calculation for the quan- 
tum system does not depend on the precise form of Hb , 
it does depend on the approximation for Hb or Ab(w). 
How can we control this dependence? There is no defi- 
nite answer to these questions, and the need to discretize 
Ab(w) is quite unsatisfactory. Our new proposition, the 
use of the CS in the CSM, avoids the discretization of 
Ab{lo). It works without a representation of Hb in the 
form ([7]) , but addresses Hb only via the spectral function 
A B (u). 

The central ingredient of CS(M) is a representation of 
Hb related to the Chebyshev expansion of Ab(u>). To 
obtain this representation, define the Chebyshev vectors 

\n) = T n {H B )S\va,c) (13) 

for n > 0. These vectors are neither normalized nor 
orthogonal. We need only the scalar products 

(0|n) = (v&c\dT n {H B ) d f |vac) = , (14) 

where fi B is the n th Chebyshev moment of the bath spec- 
tral function Ab{u>), according to the expansion (UJ). The 
Chebyshev vectors span a Hilbert space Tt c , the CS. By 
definition, H c is a subspace of the Fock space for the bath 
operators f"\ but we refer to the Chebyshev vectors 
only as abstract vectors, whose possible representation 
in terms of the is irrelevant. 

From the recurrence relation ((!]) it follows that the 
operation of Hb on TL C is given by 



H B \n) = 



|1) , n = 

(l/2)(|n-l) + |n+l» ,n^0 



(15) 



Since \ri) is a one-fermion state, d\n) is proportional to 
| vac). Together with the definition of /j, b this yields 



d\n) = dT n {H B )d)\v&c) 

= |vac)(vac|dT n (iJs)c?' l '|vac) = ii^|vac) 



(16) 



Sd\n) = d)dT n {H B )d)\ vac) 

= d t |vac)(vac|dT„(i7 B )rf t |vac) = Mn 1°) 



(17) 



The reader should note that in these equations the only 
parameters are the Chebyshev moments /j b of the spec- 
tral function Ab(w). 

Eqs. (TT4")) - (fTT|) constitute the basis of our CS approach. 
What we have achieved is that Hb is put into a form that 
makes no reference to a representation like (J7J . Still, Hb 
is given as a Hamiltonian acting on a Hilbert space. This 
is a very useful form for Hb, which can be used within 
numerically exact techniques like Lanczos or KPM and 
avoids discretization of Ab(uj). 



IV. A FERMIONIC SITE EMBEDDED IN A 
HOST 

We illustrate the use of the CS(M) with the example 
H =-AdU + H B (18) 

of an unperturbed system Hb with a perturbation 
— AdJd. This is the Hamiltonian of a single site in a 
sea of non- interacting fcrmions, like an impurity embed- 
ded in a host or lattice. Note that this example is not 
of the form ©, as it deals with a site embedded in the 
bath, while in Eq. ^ and later examples fermions leave 
the bath by hopping to the quantum system. 
Our goal is to calculate the spectral function 



A(uj) = (vac\dS(uj - H) d^vac) 



(19) 



to given Ab(oj) and A within CSM, i.e. using KPM with 
the Eqs. ([T1])-([T7|) for the CS. We can compare the results 
to the exact result 



G{z) = [Gb{z)~ 1 + A] 



(20) 



for the corresponding Green functions 

G (B ){z) = (vac\d[z - iT^-^lvac) . (21) 

To obtain the Chebyshev moments fi n of A(u>), we have 
to recursively calculate the vectors X r j(ff)tf f |vac) accord- 
ing to Eq. (fl]). The calculation proceeds in the space 
TL C , and each vector is given as a linear combination of 
Chebyshev vectors |n). With Hb according to (TlS|) and 
Ad^d according to (fl~7|) . H is represented by the matrix 



/-2A{jl b l-2A/xf 
2 
1 



-2A M f -2A//f 
1 

1 

1 



(22) 
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with H\n) = Y, m ( H )mn\m). Note that (H) mn ^ 
(m\H\n), because the Chebyshev vectors are not orthog- 
onal. Especially, (H) mn is not symmetric. Nevertheless, 
H is hermitian by definition (an explicit calculation is 
given in App. [Cjl . 

Starting with the vector |0) = GT|vac), in any step of 
the Chebyshev iteration we apply H according to (j2"2")) 
to obtain the vector T n (H)\0) from previous vectors 
T„_i(_ff)|0), T„_ 2 (ff)|0). The moment (j, n is finally ob- 
tained from the scalar product (0\T n (H)\0) according to 
Eq. (|14[) . During the iteration the index of any vector \n) 
is increased at most by one in Eq. ([22]) , and the vector 
T n (H)\0) is a linear combination of Chebyshev vectors 
|m), with m < n. The n th Chebyshev moment /i„ of 
A(lo) is therefore obtained from the first n Chebyshev 
moments /ijf of Ab(w). We have thus devised 
a computational scheme to map n moments of a given 
spectral function j4g(w) to one part Hb of H to n mo- 
ments of a spectral function A(tu) to the full Hamiltonian 
H without resorting to an explicit representation ([7]) of 
Hb in an orbital basis This is the essence of CSM. 

Note that only spectral functions Ab(u), A(cj) occur 
in the calculation. The 'missing' real part of the corre- 
sponding Green function, which is needed in Eq. (|20p . 
is implicitly accounted for by causality relations pre- 
served in the calculation. As shown in Ref. [U, we can 
obtain the Green function G(w) from moments fj, n of 
the spectral function A(uj) without invoking Kramers- 
Kronig- relations . 

A. Scaling of H and Hb 

We have so far omitted the scaling of H and Hb to the 
interval [—1, 1]. To determine the scaling of Hb we must 
choose a scaling interval Is that contains the domain of 
non-zero values of Ab{uj), i.e. Ab(uj) — for lu ^ Is (cf. 
Sec |n]). For Is = [w^ in , il> the scaling factors 

r = (^ ax - w sj/2 , s = « n + u&J/2 (23) 

give the scaling Hb = tHb + s of Hb ■ The smaller Ib , 
the higher is the resolution of j4g(w) for a given number 
of moments. Similarly, we choose a scaling interval I with 
A(lo) — Qioru^I. With scaling factors p, q determined 
from I as before, we find for the scaling of H 

H = -(H-q) = — dU+-H B + S —^- . (24) 
P P P P 

It is straightforward to introduce the scaling factors into 
Eqs. CE5]), (16]), and the matrix form ([22]) . 

Since we do not know A(lo) in advance, we have to 
rely on estimates for I. For Eq. (fT8|) . the term — Acftd 
has eigenvalues and —A. It follows that any I with 
ID JbU (Ib — A) is a possible choice. Similar estimates 
can be obtained in other situations. Alternatively we can 
determine the minimal and maximal eigenvalue of H by 
the Lanczos algorithm (cf. Sec. IVIli]) . Note that precise 



knowledge of the minimal or maximal eigenvalue of H is 
not necessary, but any bound works. It often suffices to 
obtain an estimate from operator norms. A factor of the 
order of two in the scaling is tolerable for most practi- 
cal purposes. The resolution of CSM - or KPM - still 
compares quite favorably to the resolution of other tech- 
niques. The possible loss of resolution for a too large I 
can be compensated for by increasing the number of mo- 
ments in the calculation. We show in Sec. El that we can 
even increase the resolution for A(lu) without increasing 
the number of moments for j4g(u). 

The reader should be reminded that the need to con- 
sider scaling is the price one has to pay for the high res- 
olution of the Chebyshev technique this is certainly no 
drawback of the method. Interestingly, there are situa- 
tions where the divergence of moments for wrong scaling 
can be used to our advantage: We will exploit it below 
to determine the groundstate energy of H. 

Now consider the case A = in Eq. (fig)) . We know 
that then A(lo) = Ab(w). With scaling intervals I D Ib, 
the moments fx n , /j,^ of these two identical functions are 
different since I ^ Ib- We can use CSM to calculate the 
/x n , reproducing Ab(u>) on the larger interval I (see next 
subsection) . Now assume that we choose Ib much larger 
than necessary. In principle, a possible scaling interval 
I for H might then be smaller than Ib, i.e. I C Ib- 
Will CSM then reproduce the spectral function Ab(w), 
respectively its moments, on the smaller interval II 

The answer is no: The interval I must always be larger 
than Ib, i.e. I D Ib- Then, and only then, the Cheby- 
shev iteration for A = 0, and H = Hb, is stable in the 
space Ti. c . This can be understood from the fact that Hb 
acts within H c like an operator with unity norm: The 
column sums of the matrix (f2"2"]) for A = are unity, and 
the first column has a single unity entry. Note that this 
is a statement about the norm of Hb as a matrix, and 
not with respect to the scalar product of Chebyshev vec- 
tors. If now I C j Bi then r/p > 1, and (r/p)\\H\\ > 1 
in Eq. (j2"4")) . Then, the coefficients in front of the Cheby- 
shev vectors \ri) diverge during the iteration, while the 
moments /i^ decay with the same rate. Each moment 
fi n is a sum of products of coefficients and moments ■ 
Within exact arithmetics the correct moments fi n would 
be reproduced. In finite precision arithmetics however, 
numerical round-off errors entirely ruin the output. On 
the contrary, for I D Ib, the recursion is perfectly sta- 
ble, up to any number of moments. This is an important 
point to be aware of. 



B. Numerical results 

We now show results from CSM for the model defined 
in Eq. (fig)) . We start with a semi-circular spectral func- 
tion 



A B (u;) = _V^ 2 /4-w 2 , \u\<W/2, (25) 



G 




FIG. 2: (Color online) Spectral function A(u>) for semi- 
circular Ab(w) according to Eq. (|25|) . with W = 1 and 
A = 0. The scaling intervals are I B = [— 0.6W, 0.6IF] and 
/ = [-1.2W, 1.2 W], i.e. r/p = 0.5 in Eq. $24$. The figure dis- 
plays the region around the lower band edge u; — —W/2, with 
curves for N = 2 10 ,2 13 ,2 16 Chebyshev moments, the exact 
result, and a curve with Lorentzian broadening ri/(uj 2 + rj 2 ) 
to i) = 5 x 1CP 3 . The latter curve illustrates the typical 
quality of results from the Lanczos Recursion Method. With 
N — 2 10 Chebyshev moments the exact curve is much bet- 
ter reproduced than with Lorentzian broadening. The inset 
shows the full curve, whose deviation from the exact result is 
below linewidth on this scale. The curve for N = 2 16 lies on 
top of the exact result even in magnification of the band edge, 
and demonstrates that CSM can achieve arbitrary precision. 



realized e.g. as the density of states in a Bethe lattice. W 
is the band width. The first test of CSM is the calculation 
of A{lo) for A = in Fig. H Then, A(w) = A b (lu), and 
the numerics has to reproduce Eq. (|25p on an interval / 
different from Ig. In particular, the moments \i n and 
are not identical. 

We expect strongest deviations from the numerical re- 
sult to Eq. (|25|) close to the band edges, where the spec- 
tral function behaves like a square root. For iV = 2 10 mo- 
ments, the absolute deviation |A n (Lj) — j4 e (w)| between 
numerical (A n (ui)) and exact (A e (cu)) spectral function is 
smaller than 6 x 10 -2 at the band edge, and smaller than 
5 x 10~ 5 over the remaining 90% of the band. The cu- 
mulative deviation J \A n (u>) — A e (u>)\du) is below 10~ 4 . 
Notably, the resolution is much higher than for other 
methods, e.g. using the Lanczos Recursion Method with 
Lorentzian broadening of the spectral function. Note also 
that the two-fold Chebyshev recursion - once for the mo- 
ments of j4b(w), once for A(u>) - is perfectly stable for 
any number of moments. The calculation of many thou- 
sand of Chebyshev moments is no problem, allowing for 
unprecedented resolution. 

To further test the accuracy of our approach we calcu- 
late the groundstate energy Eq(A) of an impurity in an 
infinite chain (Id), square (2d), cubic (3d) tight-binding 
lattice, and for the semi-circular density of states in 



FIG. 3: (Color online) Groundstate energy Eq(A) for an im- 
purity in a Id, 2d, 3d lattice, and with a semi-circular density 
of states. -Eo(A) is calculated with the method described 
in text, for N = 2 10 Chebyshev moments. The arrows in- 
dicate the critical A c , with A c [semi-circular] = W/4 and 
A c [3d] ~ 0.330 x W . The inset displays the error to the exact 
result for Id and semi-circular density of states, where a sim- 
ple exact result for Eo(A) is available. Already for N — 2 7 
moments we obtain results with an error below 10" 4 . For Id, 
the error decreases rapidly for larger A, since the groundstate 
is then strongly localized at the impurity, and the band edge 
singularity of the density of states less important. 

Eq. ([25| (see Fig. [3|). For sufficiently large A a bound 
state at the impurity exists, with energy outside the band 
of continuum states, i.e. Eq(A) < —W/2. In Id and 2d, 
a bound states exists for any A > 0, while in 3d and for 
the semi-circular density of states a bound state exists 
above a critical value A c . In Fig. [4] we show the impu- 
rity spectral function for A ^ 0. A(uj) has a 5-peak at 
u) = Eq if A > A c . Even if the peak is very close to the 
band edge it can be resolved within CSM by increasing 
the number of moments N. 

To find Eq(A) within CSM, we exploit the possible di- 
vergence of moments in the Chebyshev recursion. We ini- 
tially choose a scaling interval I — [cJ m i n , (jJ max ] for H so 
large that moments do not diverge. Then we vary ui m i n . 
If the moments diverge, Eq < w m i n , otherwise Eq > w m i n . 
Using e.g. a bisection algorithm on w m i n we can calculate 
Eq to high precision. The computational effort for Fig. [3] 
is independent of the lattice dimension, as only a fixed 
number of moments enter the calculation. The most de- 
manding case here is indeed Id, as the spectral function 
diverges at the band edge. Despite this particular com- 
plication, the results are extremely accurate already for 
TV = 2 10 moments. 

For our example, Eq is the groundstate energy of a sin- 
gle particle, i.e. the smallest energy with A(Eq) ^ 0. The 
weight A{Eq) can be arbitrarily small. From Eq. ([20)) . we 
see that A(uj) ^ whenever Ab{u) ^ 0. Within the nu- 
merics the domain of non-zero values of Ab(u) is fixed 
by the scaling interval Is — [^mim w max]- W e also know 
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FIG. 4: (Color online) Spectral function A(u>) for semi- 
circular Ab(lo) (Eq. [f25)). with W = 1 and finite A. The 
scaling intervals are 1b = [—(0.5 + e)W, (0.5 + e)W] and 
J = [-{0.5+e)W-A, (0.5+e)W], with a small offset e = I0" 4 . 
For A > A c = W/4, A(uj) has a pole outside the band of con- 
tinuum states (cf. Fig. [3j . The deviation of the curves with 
TV = 2 10 to the exact result is below linewidth, except for the 
pole in the curve for A = 0.4, which attains a finite width. 
The inset displays four curves for A = 0.26 and different TV. 
For this A, A(u) has a pole at E (A) sa -W/2 - 3.846 x I0~ 4 . 
Although the pole is separated from the band edge —W/2 by 
less than 10~ 3 times the bandwidth W, it can be resolved 
within CSM provided TV is large enough. It is hardly possi- 
ble to achieve a similar increase in resolution with e.g. the 
Lanczos Recursion Method. 



from the previous subsection, that moments diverge for 
^min > ^min- As a consequence, the calculation results in 



Eq < t^min- Therefore w^ in must be chosen at the lower 
band edge of Ab(cj), i.e. cj^ in = — W/2 in all examples. 
As long as \E \ > |w^ in |, especially for A > A c , no prob- 
lem occurs. Note that the situation is different at finite 
fermion density, where the groundstate energy is deter- 
mined by filling states up to a certain density n, with 
n = A{w)du). Negligible weight in A{uj) does not 
change Eq then, and the calculated Eq does not depend 
that strongly on Ig. 

From the examples shown here we conclude that CSM 
provides highest accuracy. The computational effort to 
obtain the results shown is small - a calculation for 
TV = 2 10 takes less than a second - which leaves plenty 
of room for applications to less simple problems. In the 
next paragraph we discuss how the computational effort 
can be further reduced, by adjusting the resolution of the 
two Chebyshev expansions for A{ui) and Ab(u>). 



V. PROPERTIES OF THE TRUNCATED 
CHEBYSHEV REPRESENTATION OF H B 

In the examples of the previous section we calculate TV 
moments of A(oj) from the same number of moments of 



Ab(u}). Calculations can also be performed for a different 
number of moments of Ab(ui), say M. For M > TV, the 
moments of A(u>) do not change and the calculation is 
exact, as we noted above. For M < TV however, the 
calculation is approximate and its result depends on M . 

The computational effort for given M, TV scales as 
0(MN) for time, and 0(M) for memory: Calculating 
more moments of A(uS) requires only more time, while 
using more moments of Ab(u>) requires more time and 
memory. It is therefore natural to ask whether we can ob- 
tain TV moments of A(u>) from fewer moments of Ab(u>)<, 
especially if we remember that both are calculated with 
different scaling for H and Hb- We know that the reso- 
lution of a Chebyshev expansion increases linearly with 
the number of moments, as can be deduced e.g. from the 
width (standard deviation) a = 7r/TV of the Jackson ker- 
nel, see Appendix [XJ If we demand that the resolution 
for A(u>) and Ab(w) is the same, we obtain the condition 
M /TV > r/p, with the scaling H = pH+q, H B = rH B + s 
as in Eq. As we discussed there, the factor r/p is al- 
ways smaller than unity. For the example of the Holstein 
model in Sec. IVIII A[ it is easily of the order 10 _1 . 

In Fig. Owe show A{u>) for different ratios M/N. As 
long as the condition M/N > r/p is fulfilled, we pro- 
duce an absolutely accurate A(uS) from M < TV moments 
of Ab(uj), while the computational effort is considerably 
reduced. In the worst case, for M/N <C r/p, discrete 
energy levels in the truncated bath Hamiltonian Hg are 
resolved. Even TV ^ M does not lead to erroneous re- 
sults as a negative spectral function. In particular, the 
CSM is perfectly stable for arbitrary TV and M. 

Using M < TV moments of Ab(w) in the calculation is 
equivalent to working with a truncated bath Hamiltonian 
Hg on the M-dimensional subspace Ti^f of Tt c that is 
spanned by vectors |0), . . . , \M — 1). The matrix of Hg , 
with Hjf\n) = E m {H^) mn \m), is according to Eq. jftn 
the tridiagonal M x M-matrix 



2 



(Hb )mn — — 



1 

1 

1 

1 



V 



1 

1 

1 0/ 



(26) 



We show in App.[B] that the characteristic polynomial 



of (ffjf) 

run 



det[x - {H 



B I 



2-^T M ix) 



(27) 



The M th Chebyshev polynomial Tm{x) = 
cos(M arccos x) has M distinct real roots 
for j = 1, . . . , M. It follows that 



7T(j-l/2) 



M 



(Hb )mn has M distinct real eigenvalues, hence is a 
diagonalizable M x M-matrix with real eigenvalues Xj. 
The position of the discrete energy levels resolved for 
TV M in Fig. [5] is determined by the Xj. We discuss in 
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FIG. 5: (Color online) Similar to Fig. [2] this figure shows 
the spectral function A(uj) for semi-circular Ab(u>) according 
to Eq. (f2"5")l . with W = 1 and A = 0. The scaling intervals 
are I B = [-(0.5 + e)W, (0.5 + e)W] with e = HT 4 and / = 
[-2W,2Wi, i.e. r/p = 1/4 in Eq. JM}. The calculation is 
performed for N — 2 8 Chebyshev moments of A(u>) and M 
moments of Ab(u>). As long as M/N > r/p, i.e. here M > 
2 6 , the result can not be distinguished from the (numerically 
exact) curve for M = N. For M = 2 5 , deviations occur at 
the band edges, and for M = 2 4 single spurious peaks are 
resolved. The inset displays a calculation for N = 2 11 and 
M = 2 4 . Even for this strong mismatch, a perfectly positive 
spectral function is obtained, resolving 5-peaks at uij = rxj+s 
corresponding to the roots Xj of Tm(x) (see text). 



App. [D] how the weight of the corresponding peaks can 
be obtained from the moments 



VI. LINEAR CHAIN COUPLED TO A BATH 

A different example is provided by an electron hopping 
along a finite chain of length L that is connected to a bath 
at its 'right' end (site L). The Hamiltonian 

L-l 

H = -t Y,(4+id + cjc i+ i) +H L + H B (28) 
i=i 

has the form of Eq. ([6]). Here, cj creates an electron at 
lattice site i, and the coupling between chain and bath is 
H L = -t{ScL + c\d), as in Eq. (TO]), with V = -t. That 
the bath couples to the d-orbital implies that [cf\H B ] = 
0, but [d,M,H B ] ^ 0. 

In difference to the example in Sec. IIVI (see Eq. (|T5|) ) 
the particle hops to and from the bath, hence the Hilbert 
space Ti = Hs©H c of the problem is the direct sum of the 
Hilbert space Tts of the chain and the CS He- A basis of 
Hs consists of vectors \ipi) = c-|vac), for the electron at 
site i. The vectors \i[>i) are orthogonal to the Chebyshev 



vectors |n). The operation of H is summarized in 

H\iln) = -t\ih) , 

H\i/>i) = -t|Vi+i> -A^i-i) 

for i = 2, . . . , L - 1 , (29) 
H|^> = -*|0)-#i-i), 
H\n) = -t^ L )+H B \n) , 

with the missing equations supplied by Eq. (| 1 5[) . 

The calculation of the spectral function An(ui) = 
(vac|ci5(o; — H)c\ |vac) - the spectral function to the 'left' 
end of the chain in Fig.Q]- to given bath spectral function 
A b {oj) proceeds along the lines established in Sec. IIVI In 
contrast to the example (fl~8"|) . where Eq. (fT4"|) was used, 
scalar products with the starting vector c||vac) of the 
Chebyshev iteration do not involve the fi^. 

To make the present example self-consistent we de- 
mand that Au(lj) = Ab{u:)- The self-consistent solution 
An (a;) is the spectral function of a half-infinite chain 
at its open end, i.e. the semi-circular spectral function 
from Eq. ([23]) with W — At. To obtain the self-consistent 
A\x{u)) we start with an initial guess for e.g. set- 

ting fi^ — o corresponding to A B {u) — 0, and calculate 
the moments fj, n of An(u). Note that the choice ^ = 
does not correspond to a spectral function, as the sum 
rule J A B (u>)du! = fj, B = 1 is violated. This does not 
matter for the calculation, and we could obtain the same 
effect by setting V — 0. We start a new calculation tak- 
ing the fi n just calculated as new bath moments /i^, and 
reiterate this calculation until the moments fj, n , or equiv- 
alently An(u>), are converged (see Fig- [6]) . 

We found in Sec. |W] that we must choose the scaling 
interval / for An (to) larger than I B for A B (uS), If we start 
a new iteration with the previously calculated moments 
Hn replacing /j, B we must also replace the interval I B by 
/, and consequently allow for growing scaling intervals - 
and corresponding loss of resolution - during the itera- 
tions. As an alternative, we keep the interval I B fixed 
and rescale Ah(lo) from / to I B in our implementation. 
We first construct Au(w) on the interval / from the fi n , 
then rescale it in w-space to the interval I B with a linear 
transformation uj (r/p)(uj — q) + s (cf. Eq. I24[) . and fi- 
nally feed the moments of this rescaled spectral function 
as new back into the calculation. In the rescaling we 
can check whether we throw away significant weight of 
Ah(lo) which signals a too small I B . For some exam- 
ples where we do not know the relevant interval I B in 
advance, e.g. for the Holstein polaron (Sec. IVIII Ap . we 
let the program determine a suitable I B that contains all 
but negligible weight of Ah(lu). Typically, the interval 
changes during the first few iterations and then stabilizes 
with convergence of Ah(uj). The rescaling transforma- 
tion needs two Fourier transforms (preferentially FFT) 
for the calculation of the spectral function from the mo- 
ments and vice versa, and one interpolation for the linear 
scaling. We use spline interpolation, but for not too few 
moments everything also works without interpolation. 
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FIG. 6: (Color online) Spectral function Ah(oj) for a chain 
with L = 5 sites, calculated for N = Nb = 2 10 moments. In 
the first iteration, An(ui) consists of L peaks. The number 
of peaks is increased by L per iteration (see curve for 2 itera- 
tions). Peaks merge if their distance is smaller than the reso- 
lution (see curve for 17 iterations), until Ai\(u) is converged 
to the self-consistent semi-circular solution after a sufficient 
number of iterations (here about 30 iterations). The scaling 
intervals are I B = [-(W/2+e), (W/2 + e)], with e = 10~ 4 , and 
I = [-W,W], with W = it. Note that the number of itera- 
tions until convergence depends on the number of moments, 
i.e. the resolution. For eg. N = 2 7 moments, convergence is 
obtained after 6 iterations. 



The reader should note that the rescaling transfor- 
mation - which shrinks the w-interval of the Cheby- 
shev expansion - is a linear transformation of moments 
{l*n} !— > A calculated fi^ depends on every 

|U„. In contrast, the CSM of Sec. IIVI implements for 
A = in Eq. (|18p a transformation {//^} i— > {fJ, n }, that 
blows up the w-interval (cf. Fig. An inspection of 
Eqs. (|14p. (|15p shows that this transformation is also lin- 
ear. For A ^ 0, it ceases to be linear as the fjg also oc- 
cur in H itself (cf. Eq. by virtue of Eq. (TTT]). The 
transformation has the property that each calculated [i n 
depends only on //^ with m < n. We now understand 
a second reason why the Chebyshev iteration has to be 
unstable for Ib <f. I- If the Chebyshev iteration could be 
used to shrink the w-interval - that is for I c Ib ~ every 
calculated moment (p n ) had to depend on every moment 
supplied to the calculation (/x„ ), in contradiction to the 
properties of the recursion. 

We discussed in Sec. IVlhow to calculate with CSM N 
moments /i„ from M moments /i^ for M < N, depending 
on the scaling intervals I and Ib- Do we waste part of 
the numerical results if again only M < N moments /z^ 
are calculated in the rescaling transformation? Obviously 
not: For a smaller interval less moments are needed to 
obtain the same resolution. In both directions, via CSM 
or the rescaling transformation, the number of moments 
should be related by the estimate M/N ~ r/p. 



VII. CONCURRENT DYNAMICS 

The general Hamiltonian Eq. ([6]) operates on a product 
space H = Hs ®Hb. H$ {Hb) operates on Hs (Wb), 
and Hl links the two spaces. In CSM, Hb = H c - We can 
generally write the Hilbert space of H in such a way, with 
Hs and Hb as Fock spaces, but additional constraints 
may single out a subspace. In the example of the linear 
chain with the restriction to a single electron, H is the 
direct sum of Hs and H c - 

Assume that system and bath do not couple, i.e. Hl = 
and H = Hs + Hb- As Hs and Hb commute, the 
spectral function A(uo) = (ip\6(u> — H)\ip) to a product 
state \ip) = \ips) ® IV'b) is the convolution of the spectral 
functions of Hs and Hb ■ Explicitly, 

A(Lu)=J2\(^s)\ 2 A B (uj-e a ) , (30) 

a. 

where the sum is over eigenstates \a) with eigenvalue e a 
of Hs, and Ab{u) = (i/jb\S(u) — Hb)\iPb)- To evaluate 
this expression we must diagonalize Hs in advance, to 
obtain its eigenvalues and eigenstates. It is therefore dif- 
ficult to use Eq. (f3"0"]) for even moderately complicated 
Hs- Within CSM we calculate A(u>) without diagonal- 
ization of Hs- The convolution in Eq. (|30p is implicitly 
performed in course of the Chebyshev iteration. This fea- 
ture is essential for all situations where system and bath 
degrees of freedom evolve in parallel. 

We want to illustrate this point for the Hamiltonian 

H s = -V^o( &t +6)c t c + w & t 6 (31) 

of a bosonic site, the independent boson modet^. It is 
a simple model for electron-phonon coupling of localized 
electrons, if the bosons parametrize the elongation 
of an ion that produces an electric field which shifts the 
energy of an electron at the site (c'). The groundstate 
of the bosonic part of Hs in presence of a fermion, to 
energy Eq = — e p , is the coherent state 




For Ab(u) we assume a semi-circular spectral function 
as in Eq. (|25|) . and prepare the system in the state \ip) = 
| cob.) (g> c^vac). We calculate the spectral function 



A(w) = {ip\c)d8{uj -H + E )Sc\^) . (33) 

Physically, this corresponds to a sudden excitation of the 
electron from the bosonic site (c-orbital) to a continuum 
of states given by Hb, which is related to X-ray absorp- 
tion of a localized electron^. We know that we should 
obtain A(ui) as a sum of semi-circular spectral functions 
shifted by multiples of luq, weighted with the Poissonian 
distribution of bosons in |coh). The lowest band is cen- 
tered at e p , which is the energy to remove the electron 
from the bosonic site. In Fig. [7] we show A(ui) calculated 
by CSM. The numerical result perfectly agrees with the 
expected outcome. 
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7: Spectral function A(ui) from Eq. (|33[) , f° r = 1, 
€ P = 4 in Eq. fl3T|, and VF = 0.5 in Eq. J25J). We use N = 
2 13 moments and M = 2 s bath moments. For maximally 
n& = 25 bosons in the calculation, the scaling of H and Hb is 
r/p < W/(n b o;o) = 0.02. According to Sec. El it is sufficient 
to use about M = 160 bath moments. The inset displays 
in magnification, how accurately the semi-circular spectral 
function is resolved for a single subband. 



VIII. A BOSONIC SITE COUPLED TO A BATH 
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FIG. 8: (Color online) Average number of bosons (6^6), oc- 
cupation probability of the impurity site (e'e), and ground- 
state energy Eo in dependence on e p , for the model (|34|l with 
A/W = — 1 and a semi-circular Ab(uj) with W = 1 according 
to Eq. (|25[) . Calculations have been performed with Lanc- 



zos/CS for M = 2 9 bath moments. Already M = 2 mo- 
ments produce accurate results away from the phase transi- 
tion. Close to the phase transition, we can increase the accu- 
racy easily by increasing M. Note that M does not need to 
exceed the number of Lanczos iterations. 



We assumed in the previous section that Hl = 0, so 
that system and bath degrees of freedom do not mix. 
Nothing changes to the applicability of our approach if 
this restriction is abandoned, as the next examples show. 
In contrast, the calculation of such 'mixed' dynamics is 
the intended application of CSM. 

Let us combine the example from the previous section 
with the example from Sec. HVl We get the Hamiltonian 

H = - Ac f c - v /e^Tj u '(& t + 6)c f c + u tfb 

-t{Jc + cU) + H B '" 4l 

of a bosonic impurity (c-orbital) coupled to a fermionic 
bath, e.g. a lattice. Note that the impurity site is coupled 
to the bath via a term Hl = —t(<vc + c'd) as in Eq. ([5]) 
or Eq. ((28)) . We could also study a model with an bosonic 
impurity embedded in a host similar to Eq. (|18[) . but the 
present form is convenient for the study of the Holstein 
model (Sec. I VI II AD. 

For A < 0, the impurity is repulsive and acts as a static 
barrier for electron motion. Two competing mechanism 
determine the groundstate: On the one hand, the energy 
of the impurity state is increased by —A, favoring a de- 
localized groundstate. On the other hand, the formation 
of a localized polaron at the bosonic impurity lowers the 
energy of the electron roughly by e p . A localized impu- 
rity state occurs if the loss in kinetic energy is overcome 
by the gain in potential energy. We know from Sec. IIVI 
that this happens - for an attractive impurity - if A is 
larger than a critical A c . By a rough estimate, we expect 
here a localized groundstate for e p > A c — A. 



To address this issue we must calculate the ground- 
state of the model (|34|) . In Sec. II V Bl we explained how 
to determined the groundstate energy Eq using Cheby- 
shev expansions, testing for the divergence of Chebyshev 
moments for different scaling of H. This time we must 
obtain the groundstate itself, not just its energy, to be 
able to calculate correlation functions for the degrees of 
freedom of the bosonic site. To calculate the groundstate 
we use the Lanczos algorithm. This requires our ability 
to perform two operations: To apply H to a vector, and 
to calculate the scalar product between vectors. For Hb 
defined on the CS Tt c , Eq. (fT3)) defines the application of 
H B to a vector, and Eq. (or Eq. JC3]) in App. EJ) 
gives the scalar product. Note that for a successful ap- 
plication of Lanczos the moments /i^ should be modified 
by attenuation factors (Eq. (|A5|I ), cf. App. lAlandlO 

In Fig. [5] we show the groundstate energy Eq, the oc- 
cupation probability (e'e) and the average number of 
bosons (b^b) calculated with the Lanczos algorithm and 
the CS. As in Sec. IIV Bl the computational effort is in- 
dependent of the dimension. We used a semi-circular 
Ab{oj) here, for a 3d lattice the results are qualitatively 
the same. While Eq is a smooth function of e p , (e'e) and 
(b'b) signal the phase transition from a delocalizcd elec- 
tron to a localized polaron at a critical ep\ For luq — > oo, 
£p converges to the value A c — A of the simple estimate. 
For luq — * 0, the phase transition becomes more pro- 
nounced, as a precursor of the first-order transition in 
the adiabatic limit ujq — 0. Note that in contrast to 
the (Holstein) polaron problem, this model has a phase 
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transition also for finite loq. 

It is perhaps surprising that the Lanczos algorithm, 
which constructs an orthonormal basis, works in combi- 
nation with the CS construction, which is based on non- 
orthogonal vectors. For H = Hb the Lanczos algorithm 
constructs an orthonormal basis of H. c and is equivalent 
to Gram-Schmidt orthonormalization of the Chebyshev 
vectors \n). It is known that the Gram-Schmidt proce- 
dure is prone to instabilities, namely loss of orthogonal- 
ity. The same is true for the Lanczos algorithm, which is 
a problem for the calculation of spectral functions. The 
calculation of the groundstate however does not require 
to construct an orthonormal basis of H c , and loss of or- 
thogonality is not a severe problem. For the model (fTS)) 
or (|34[) we can force the Lanczos algorithm to fail, if 
A = e p = and w^ in <C — W/2 for the lower bound 
o;^ in of Ib- For these parameters, the groundstate is a 
linear combination of Chebyshev vectors |n) with energy 
Eq > 0J^ in , and the Lanczos algorithm will not converge. 
Remember also the discussion in Sec. IIVBI concerning 
the calculation of Eq(A) for a single fermion, where we 
noted that we must guarantee Eq < w^ in . In practice, we 
never encountered a problem with the Lanczos algorithm 
in combination with the CS construction. For Fig. ([5]), 
Lanczos converges fastest for large e p , when the ground- 
state is localized at the impurity site. 

A. Holstein model 

The solution of the model is related to the self- 
consistent solution of the Holstein model within DMFT. 
The Holstein models is a standard model for electron- 
phonon-coupling. Its Hamiltonian 

H = - *w °i c i ~ y/tpUo ^2(b\ + bi)c\ Cl + ujq b\bi 

lj i i 

(35) 

contains a local coupling of the electron density (cjcj) 
to dispersionless optical phonons (&J), in addition to the 
kinetic energy of electrons modelled by the hopping term, 
and the kinetic energy of phonons with frequency luq. 
The polaron shift e p is the groundstate energy of the 
model for t tj = (cf. Sec. ED- 

The solution of the Holstein model, especially for spec- 
tral properties, is still a demanding problem (for a recent 
review see Ref. [IB). We successfully used Chebyshev ex- 
pansions to obtain spectral functions or the optical con- 
ductivity for finite system s 17 ' 18 . Here we use DMFT to 
relate the Holstein model to the model (13~3| where a cou- 
pling to a bath occurs. 

Within DMFT, the (fc-integrated) spectral function 

A(u>) = (vac|c 4 <S(w - ff)4|vac) (36) 

is obtained as the self-consistent solution of an impurity 
model with a single interacting site^ For the Holstein 
model p5p . the impurity model is just the model of a 
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FIG. 9: (Color online) Spectral function A(u) and imagi- 
nary part of the selfenergy for the Holstein model (|35[) within 
DMFT, for a semi-circular non-interacting density of states 
with W — 1 (Eq. (|25p ). Calculations were performed for 
M = 2 11 bath moments. The phonon frequency is luq/W = 

0. 25, the coupling strength e p /W = 0.25 (upper panel) and 
e p /W — 0.75 (lower panel). Note the different scale for 
ImE(oj) in the lower panel, where a pole of £(oj) separates 
the lowest band. The inset displays a magnification of the 
lowest polaron band and the pole. 

single bosonic site for A = 0. If we assume a semi- 
circular density of states for the non-interacting problem, 

1. e. A(u) has the functional dependence of Eq. (|25[) for 
e p = 0, self-consistency is established by A(u) — Ab{w). 
For a different non-interacting density of states this re- 
lation is more complicated, but the difference is not rel- 
evant in our context. We have shown in Sec. IVII how 
a self-consistent calculation is performed. The compli- 
cated part of the DMFT calculation, to obtain A(u>) for 
the Hamiltonian ([34)) to given Ab{u>), is solved by an 
application of CSM. 

In Fig. [5] we show A{uS) for two sets of parameters. For 
e P /W — 0.75, a polaron has formed as a new quasipar- 
ticle, which results in several separated bands in A(u>). 
We also show the imaginary part of the DMFT selfenergy 
£(c<j). In our example, it can be obtained in the form of 
a Green function 

E(cj) = (vac|c6(w - i2i) _ Vc + |vac) , (37) 

where Hi = (1 — P$)H{1 — Pq) is projected onto the 
subspace orthogonal to the boson vacuum |vac), i.e. 
Po = | vac) (vac | . The projection guarantees that only 
irreducible diagrams contribute to S(w). We calculate 
£(u;) once A(u>) is converged. When a polaron has 
formed, the polaron bands are separated by a pole in 
Similar to Fig. [U a pole close to a band has to 
be resolved, which requires high resolution. Note that 
ImE(u) is zero for the lowest polaron band, where emis- 
sion of a virtual phonon is energetically forbidden. 

In Fig.[TO]we show a calculation for a Id chain. Since 
DMFT is constructed in the limit of high dimension, only 
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k=n/2 



FIG. 10: (Color online) Spectral function A(u>) for the Hol- 
stein model on a Id chain (W — At with t = 1) within 
DMFT-approximation, for lu /W = 0.25, e p /W = 0.5. The 
lower panel shows the fc-resolved spectral function A(k,u) — 
A°(k,u) - £(w)), where A°(fc,a>) = S(co - e fe ) with e fc = 
— 2tcosfc is the spectral function of a Id tight-binding chain, 
and E(w) the DMFT selfenergy (Eq. 157). Calculations were 
performed for M = 2 1 bath moments. 



the basic features of polaron formation in Id are correctly 
described. With CSM, this calculation could be extended 
to a cluster of bosonic sites, using one of the recently 
developed cluster extensions 19 of DMFT. 



B. Comparison to the analytical solution 



For the Hamiltonian l(34jl. hence for the DMFT- 
solution of the Holstein model (|35p . an explicit solution 
for the spectral function A(ui) can be obtained as a con- 
tinued fractionSSiSI (CF) 



A(lo) 



u - t 2 A B (uj) 



e p uo 



lu - t 2 A B (uj - u ) 



(38) 

We can construct the CF from a formal Lanczos re- 
cursion. The recursion starts with the state cHvac) and 



consecutively produces the states (tf / Vntyrf \vac) , which 
form an orthonormal basis of the Hilbert space Tls of 
the bosonic impurity site. In this basis, the matrix of H 
(Eq. (|M|) ) is tridiagonal, and A(lj), which is the (1,1)- 
element of this matrix, can be expressed as a CF. 

To treat the coupling to the bath in Eq. ([5~4")) , we use 
Eq. ([50)) . If the electron is in the bath, the states of the 
bosonic site evolve by the Hamiltonian H' s = u^b. In 
the n th level of the CF ([55)) . which corresponds to the 
excitation of n bosons, we must therefore insert Ab{uj) 
with energy shift Ab{lo — ruvo). 



For this example, the (Lanczos) recursion leading to 
the CF creates every eigenstate of uJi^rb one by one. 
Otherwise, if eigenstates were mixed during the recur- 
sion, linear combinations WnAsi^J — nu>o) would oc- 
cur instead of Ab{u — nwo)- The weight w n had to be 
determined during the recursion, and depended on the 
parameters. This is one of the many obstructions that 
prohibit a generalization of Eq. ([55)1 to other models, e.g. 
with a cluster of interacting (bosonic) sites or different 
electron-phonon coupling. Within CSM these consider- 
ation are pointless, as its application does not rely on 
our ability to solve one part of the Hamiltonian prior 
to the actual calculation. Note also, that the general 
CS construction allows to calculate the groundstate and 
correlation functions that cannot be expressed in a form 
like O- 



IX. TIME EVOLUTION OF A WAVE PACKET 
ON A CHAIN COUPLED TO LEADS 

In the preceding sections we used the CS construction 
for the calculation of spectral properties within CSM. It is 
an essential advantage of Chebyshev techniques, and the 
CS construction, that they can be easily adapted to new 
problems. In this section we treat the time evolution of a 
wave packet on a chain coupled to leads. To describe the 
coupling to a lead we can use the CS construction from 
Sec. IIIII without change. 

Within Chebyshev techniques, time evolution of a vec- 
tor \ip(t)) = e~ lHt \tp(0)) is determined from the Cheby- 
shev expansion of the time evolution operator— 



-iHt 



Co 



(39) 



The expansion coefficients are given by Bessel functions 
- 1 T n (x)e- ixt 



-i 7rVT 



= (-i) fl J B (t) , 



(40) 



where J n (t) denotes the Bessel function of order n. As 
usually, we omit the scaling of H. Since J n {t) decays 
rapidly for n 3> t, the infinite series can be truncated to 
obtain \ip(t)) with high precision. An adequate choice for 
the number of moments is given by N > 1.5t. 

We apply the Chebyshev time evolution to the propa- 
gation of an electron along a chain of length L, which is 
coupled to a bath at site L. The chain geometry is iden- 
tical to the example studied in Sec. lVIl and the Hamilto- 
nian is given by Eq. {28J. The calculation of T n (H) |^(0)) 
proceeds along the lines established in previous sections, 
and \ij){t)) is obtained from Eqs. l[39"]). ([40)1 . 

At t = 0, we place a Gaussian wave packet of width a 
and momentum K centered at site mo = L/2, 



Mo)) 



-(m-rao) /2cr c f 



vac) 



(41) 
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FIG. 11: (Color online) Time evolution of a Gaussian wave 
packet Eq. (gTJ), witn 2cr 2 = 25 2 and K = 1. The time step 
between two curves is St = 40. The group velocity of the wave 
packet is given by the dispersion — — 2tho P cos k as v = 
(<9efc / ! dk)\n ■ Per time step, the wave packet moves v x St ~ 67 
lattice sites. The three panels correspond to the situations 
described in the text. 



(we omit normalization here), and let it evolve in time 
(see Fig. QT]). Without a bath, the particle is reflected 
at the right end of the chain and returns, moving to the 
left. 

We now add a bath to the right end, whose spectral 
function is given by Eq. (f25|) , with W — 4t. Since 

this is the spectral function of a half-infinite chain at its 
open end, also the full Hamiltonian H describes a half- 
infinite chain. Physically, the open boundary at site L is 
removed by coupling to an infinite lead. In contrast to 
the previous situation, the particle will not be reflected 
at site L, but propagates into the bath (or lead). As 
the middle panel in Fig. QT] shows, no spurious reflection 
occurs. 

In the Hamiltonian H, sites 1 . . . L are explicitly in- 
cluded while sites > L + 1 of the lead are realized through 
the bath. In this way, we can use our Chebyshev ap- 
proach to realize 'transparent boundary conditions ' 23 ! 24 
that mimic an infinite system with a finite number of 
lattice sites explicitly treated. Transparent boundary 
conditions are implemented in Rcfs. I23ll2~4l by a modi- 
fication of time propagation algorithms like e.g. Crank- 
Nicholson. Within the CS construction, modified bound- 
ary conditions are implemented independently of the ac- 
tual calculation. For time propagation we use the same 
Eqs. (fl"4]) - (fl"6|) as for the calculation of spectral proper- 
ties. This permits us to use a variety of algorithms, like 
we used Lanczos algorithm in Sec. IVIIII We could also 
use the CS in combination with Crank-Nicholson, but 
for the time-independent Hamiltonians considered here 
Chebyshev time propagation is most efficient. 

Transparent boundary conditions correspond to a par- 
ticular choice for Ab(uj). We have the freedom to choose 
any Ab(u>)- For example, if Ab(oj) is the spectral func- 



FIG. 12: (Color online) Time evolution of a Gaussian wave 
packet with same parameters as in Fig. 1111 Time is shown on 
the vertical axis. A bath which mimics a chain segment of 500 
sites is placed between sites 500 and 501 of the chain (vertical 
dashed line). The wave enter the bath and reemerges after a 
delay At = 500/n ~ 300 (upper panel). After transmission 
through the bath the wave packet is 500 sites behind the wave 
packet on a simple chain, whose propagation is shown in the 
lower panel for comparison. 

tion of a 300 site chain at its end, and L = 700, H de- 
scribes a chain of length L + 300 = 1000. As before, 
the particle is reflected and returns moving to the left 
(see bottom panel in Fig. [IT] in comparison to the top 
panel) . The point is that the reflecting barrier is now re- 
alized through a bath with carefully chosen and 
not by lattice sites 701 ... L explicitly included in the 
Hamiltonian. We might call these 'reflecting boundary 
conditions'. 

For Fig. [T3 we slice the chain in two parts which are 
connected by a bath. Similar to the previous example, 
the bath shall replace a finite segment of the chain. Since 
the bath has now two entry points - at its junction to 
the left and to the right part of the chain -, this example 
involves off-diagonal bath spectral functions. Let dy 2 
denote the operators to the two entry points. Following 
Sec. lIIIl we need two different types of Chebyshev vectors 
\nb) = T n {H B )d[\va,c) for b = 1,2. Eq. (fTTJ)) generalizes 
to dalrib) = /x^' a6 |vac), where the four different types of 
moments correspond to the four different spectral 

function characterizing the bath: The diagonal Af^u), 
A^us) for moving into the bath and back, and the off- 
diagonal Af 2 (LL>), A 2 i{iS) for moving through the bath. 
The latter ones give the energy-dependent transmission 
rate through the bath. Only for non-zero Af 2 (ui) , A 2X {lo) , 
the bath connects the two parts of the chain. In Fig. [12] 
the spectral functions of the bath are these of a finite 
chain segment. We see how the wave moves into the 
bath at its left side, passes through and reemerges at the 
right site. The bath perfectly mimics the transmission 
through a chain segment. As we noted above, we have 
the freedom to change the behavior by a different choice 
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of the bath spectral functions A^ b (uo). 



X. SUMMARY AND OUTLOOK 

In this article we introduced the Chebyshev space 
method (CSM) for the treatment of degrees of freedom 
with non-trivial dynamics. In our examples the degree 
of freedom is given by the operator <v which creates a 
fermion in a bath, and the dynamics of d) is specified 
by a spectral function We demonstrated for var- 

ious examples, how CSM yields extremely accurate re- 
sults with modest computational effort. For the example 
in Sees. IIV Bl IVIIII the computational effort is even inde- 
pendent of the geometry or dimension. 

A particular advantage of the Chebyshev space (CS) 
construction is that it still provides a Hamiltonian Hg 
acting on a Hilbert space, albeit it is an abstract space 
without a direct physical counterpart. We can use this 
form of Hb in different situations, like in the Lanczos 
algorithm or for time propagation. We have therefore 
obtained a numerically exact Hamiltonian treatment of 
degrees of freedom with non-trivial dynamics, which docs 
not rely on a discretization of spectral functions. To the 
best of our knowledge, this has been achieved for the first 
time. 

The concept of a CS and an ancillary Hamiltonian Hb 
acting on this space is very general. It makes no as- 
sumption about the Hamiltonian H$ for the quantum 
system, and can be extended to various problems men- 
tioned in the introduction. We have obtained preliminary 
results for the Holstein model within cluster extensions 
of DMFT, but the issue needs further exploration from 
the physical point of view. 

An especially important topic is the extension to in- 
teracting fermions at finite density, or to bosonic baths. 
With a single fermion, the bath is in the vacuum state 
| vac) if the fermion is removed by application of the op- 
erator d. This simplification is no longer true for finite 
fermion density, when adding and removing a fermion 
creates particle-hole pairs in the bath, which initially was 
prepared as the particle-hole vacuum (or Fermi sea). We 
have derived a formalism that deals with particle-hole 
pairs in the context of CSM. The CS construction is used 
in this formalism without change. A similar formalism for 
bosonic baths can be derived. As a first result, this for- 
malism allows for the calculation of the core level spectral 
function in the X-ray absorption problem*^. The spec- 
tral function can be calculated exactly in this case. It 
is equivalent to the Anderson model with one immobile 
spin species, but the extension to the full Kondo problem 
is much harder. The presentation of this formalism is left 
for a future publication. 

The CS construction can be also used in the context 
of (diagrammatic) Green function techniques. Any Feyn- 
man diagram contains energy-dependent Green functions 
which represent a degree of freedom with non-trivial dy- 
namics. The difficult task is to sum up a huge num- 



ber of Feynman diagrams. Note that we can interpret 
the self-consistent calculation of the Holstein model in 
Sec. IVIII Al as the exact summation of a certain class 
of Feynman diagrams. The bosonic impurity model de- 
fines a set of skeleton diagrams for this problem, and 
imposing self-consistency corresponds to replacing bare 
Green functions ('thin lines' in a diagram) by renor- 
malized Green functions ('thick lines'). In this exam- 
ple, diagrams are selected by a geometric rule, but with 
some modifications a different selection rule can be im- 
plemented. The CSM provides the link between exact 
numerical techniques for Hamilton operators, and gen- 
eral approximation schemes for Green functions. 

In conclusion, we believe that the CSM introduced in 
this work is a powerful addition to existing numerical 
techniques in theoretical physics or chemistry. It sub- 
stantially enlarges the field of applications of Chebyshev 
techniques and keeps their advantages. The results we 
obtained in this work are promising for successful appli- 
cations to more complicated problems, and the possible 
combination of CSM with other techniques shall prove 
fruitful. The further development of CSM and its appli- 
cation to the study of physical problems mentioned in 
the introduction is the subject of current research. 

Acknowledgment. We are indebted to G. Wellein 
for helpful comments during the preparation of this 
manuscript. 



APPENDIX A: CHEBYSHEV EXPANSIONS AND 
THE KERNEL POLYNOMIAL 

In this article we repeatedly approximate a function 
/ : [-1,1] -> 1 by a finite Chebyshev series 



In{x) 



N-l 



MO 



2 ^2 ^nT n (x) 



where the moments /i n are given by 



Mn = J f{x)T n {x) dx 



(Al) 



(A2) 



A central question is how good f(x) is approximated by 
/jv- Functional analysis teaches us that the /at converge 
to / in the integral norm 



\\f{x)-f N {x)\\ = 



{f{x)-f N (x)f 

7Tvl — X 2 



dx 



1/2 



(A3) 



corresponding to the scalar product used in ([2]) , provided 
||/(x)|| is finite. For practical purposes, this convergence 
property is too weak, and we demand uniform or at 
least pointwise convergence. We do not try to resolve 
the issue of convergence here, but state that for suffi- 
ciently smooth f(x) the series will converge uniformly 
to / on any closed sub- interval of [— 1 , 1] that excludes 
the end-points ±1. Convergence at the end-points is of 
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course spoiled by the divergence of the weighting function 

In physical applications, we want to expand spectral 
functions which are not necessarily continuous, and per- 
haps contain some <5-peaks due to (quasi-) particle states 
with infinite lifetime. In this case, the series cannot con- 
verge uniformly - otherwise its limit / had to be con- 
tinuous itself -, but the situation is even worse. The 
infamous Gibbs phenomenon ruins convergence at all: In 
plain words, the series (|A1[) fails to converge in the vicin- 
ity of a discontinuity, e.g. a jump of the function, instead 
shows rapid oscillations whose height does not decrease 
as the number of terms N — > oo. The Gibbs phenomenon 
is a severe obstacle for practical applications. 

Fortunately, the problem arising from the Gibbs phe- 
nomenon is solved for Chebyshev expansion a ' ' ' ■ It is 
possible to specify, prior to the calculation of the fi n , a 
set of attenuation factors <?„ , n = 0, . . . , N — 1 for every 
N, such that the modified approximants 



Sn = 



1 



ttVI 



N-l 



N 

3o fJ-o 



2j29n»nT n (x) (A4) 



n=l 



do not show the Gibbs phenomenon but approximate a 
wide class of functions in a favourable way. Plainly spo- 
ken, the attenuation factors damp out high frequency 
oscillations that otherwise lead to spurious results in the 
finite series. Instead of a more precise formulation we 
show an example in Fig. 1131 The best choice in many 
cases are factors g^ derived from the so-called Jackson 
kernel^, leading to the explicit expression 
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N 



(N 



1) cos 



7rn 
N+l 



N+l 



N+l 



N+l 



(A5) 



Modifying moments with these attenuation factors is re- 
lated to convolution with an almost Gaussian peak of 
width a — n/N (cf. Ref. 5). An important property of 
the modification with these attenuation factors is that 
the modified approximants /at are positive whenever the 
function / is. This is not true for the unmodified approx- 
imants (cf. Fig. [13]). Besides Eq. (|A5|) . other suggestions 
for attenuation factors exist, see Ref. i for a different 
choice in the context of Green functions, with an appli- 
cation in Ref. [2f| 

Note that the factors g„ in Eq. (|A5|I have the property 
that g„ — > 1 for N — > oo. We can therefore always use 
the attenuation factors g^ without losing information. 
For functions with discontinuities multiplication with at- 
tenuation factors enforces the decay of moments relevant 
for good approximation properties. For smooth f(x) the 
moments itself decay rapidly, and multiplication with the 
attenuation factors does not change the result for not too 
small N. 

The introduction of attenuation factors is the essence 
of the kernel polynomial method (KPM). Within KPM 
resolution for spectral functions is already obtained with 
a fairly small number of moments, thus allowing for ac- 
curate calculations with moderate demands on computa- 
tional time or memory. The possible resolution is much 




FIG. 13: (Color online) Shown is the spectral function to the 
Hamiltonian (|18[) and A = 0.4. The solid curve is calcu- 
lated with attenuation factors from Eq. (I A5|) and is identical 
to the curve in Fig. [4] The dashed curve is calculated with- 
out attenuation factors. The inset displays the curve from 
Fig.Uto A = 0.26 and N = 2 10 , again with (solid) and with- 
out (dashed) attenuation factors. These two curves illustrate 
the remark from the text: After multiplication with attenu- 
ation factors the reconstructed spectral function seems to be 
smooth. But the unmodified moments correspond to a dis- 
continuous spectral function with a pole close to ui — —0.5. 
The slow decay of moments of this spectral function results 
in oscillations. The pole is resolved only by increasing N and 
keeping the attenuation factors (cf. the curves in Fig. [4] up to 
N = 2 16 ). 



better than known e.g. from the Lanczos algorithm 
where peaks are commonly broadened by a Lorentzian 
function. Furthermore, the resolution of the Chebyshev 
expansion is uniform over the energy interval, and does 
not deteriorate towards the center of the spectrum. 



APPENDIX B: DETERMINANTS AND 
TWO-TERM RECURRENCES 

We formulate a well-known determinant identity re- 
lated to two-term recurrences. Define, for each n and 
given cij , bj ,Cj, a tridiagonal n x n- matrix 



/ Oi c\ 
bi a 2 c 2 

b 2 a-3 



\ 



On 

bn 



\ 



-2 C n -2 
-2 O n _l C„_l 
b n -l 0<n ) 



(Bl) 



Expanding in the last row yields 

det A n = a n det - 6„_ic„_i det A„_ 2 • 



(B2) 
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From this relation it follows that the characteristic 
polynomials Xn — det (a; — A n ) obey the two-term re- 
currence 



Xo = 1 , Xi = x ~ a i j 

Xn+l ={X~ a n+1 )Xn ~ KCnXn-l ■ 



(B3) 



Conversely, for polynomials P n defined by a two-term 
recurrence 



P = l, Pl=Oil, P, 
this result implies that 

/a±x d\ 

P.„ = det 



71-2 



(B4) 



C\ a 2 x d 2 

c 2 a 3 x d 3 
C3 04 a; 



V 



(B5) 



where c„, d n is chosen in such a way that c n d n = b n , e.g. 
for positive b n we may choose c„ = d n = \fb~^. By multi- 
plying the j th column - or equivalently, j th row - by 1 /aj , 
P n can be expressed as the characteristic polynomial of 
a tridiagonal matrix times a prefactor a j- 
In particular, the Chebyshev polynomials fulfill 



T n (x) = det 



/x 1 
1 2x 1 
12x1 
1 2x 

V 



(B6) 



Multiplying the second to the last column by 1/2, we find 
the result T M {x) = 2 M ~ 1 det(x - Hg) used in the main 
text. 



APPENDIX C: HERMITICITY OF H B 

We first prove that Hb is hermitian by calculating ex- 
plicitly that Hb is equal to its adjoint. The product of 
two Chebyshev polynomials satisfies the identity^ 



T m (x)T n (x) = (T m+n (x) +T m _„(x))/2 



(CI) 



For notational convenience, we define here T^ n (x) — 
T n {x). It follows for Chebyshev vectors |n) = T„(iJg)|0) 
that 



(m\n) = (/x„+„ + Mm-„)/2 , 



(C2) 



where we again define /i B n = /i B . The scalar product of 
two vectors \a) = J2n=o a n\n), \b) = J2n=o § iven as 
linear combinations of Chebyshev vectors is now found 
as 



Since fx B _ n - ■■ L] 



P'n-mi W ^ n the notational convention 



adopted above, we have (a\b) = (b\a)*, as it must be 
(note that the moments are real). 

To calculate a scalar product with Hb we use that 



1 



H B \a) = -J2an(\n-l) + \n + l)), (C4) 
n=0 

where the summand for n = is ao( \ — 1) + |l))/2 = ao|l) 
in accordance with Eq. (|26p. Inserting this equation we 
find 



(b\H B \a) 



2 b* m a n ((m\n-l) + (m\n + l)) 



4 ^ 

rn.n—Q 



(C5) 

i * ( B 1 B v ' 



, B 1 B \ 

"r Pm+n+l "+" H'm-n-l) • 

We can condense this equation to matrix-vector form 
1 °° 

(b\H B \a) = - ^ Ki a nM mn , 

m,n—0 

= Mm+ri— 1 Mm-n+l Mm+n+1 Mm— n— 1 ■ 

(C6) 

Since M m „ depends only on the sum and difference of 
to, n, M m „ = M„ m = M* m (the moments are real, and 
= Mf )• I* follows that (b\H B \a) = (a\H B \b)*, and 
Hb is hermitian. 

We derived the hermiticity of Hb without imposing 
any constraint on the moments. This may be a bit 
puzzling, because the moments should correspond to a 
non-negative spectral function Ab{oj) if Hb is hermi- 
tian. The puzzle is resolved by observing that the scalar 
product (|C3|) has to be positive (semi-) definite. This re- 
quirement imposes constraints on the moments that are 
equivalent to a non-negative Ab{oj). To mention one, 
(n\n) = (fi B + fi B n )/2, so it must fi B n > —fi B . Consider 
for example Hb = 0, equivalent to Ag(w) = S(oj). We 
then have fx B n = (—1)", H B n +i — 0, i.e. equality holds 
in this inequality. The scalar product (|C3[) is positive 
semi-definite, but not positive definite, as |n) = and 
(n\n) = for odd n. 

We can characterize the properties of the scalar prod- 
uct (|C3|) better, if we consider Chebyshev expansions 



f( x ) = fm.T m (x) , g(x) = g m T m (x) (C7) 



m=0 



m— 



of functions /, g : [— 1, 1] — ► M, and the associated vectors 



I/) 



f m \m) , \g) 



g m \m) 



(C8) 



(%)= 51 fo r>n(M™+„+M™-J/2 



(C3) 



in H c - With Eqs. (|CTj) . © we find for A B (co) as in 
Eq. Q f^T^T^)^)^ = ( M * +n + J/2, 
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which implies 



f(x)g(x)A B (x)dx=(f\g) 



(C9) 



The scalar product (|C3[) in 7i c therefore corresponds to 
a scalar product of functions [—1,1] — > R, given as the 
integral on the left hand side of this equation. We con- 
clude, that for continuous Ab(u>) the scalar product (|C3|) 
is positive semi-definite (positive definite) if and only if 
Ab(u>) > (and Ab(u>) does not vanish on an open in- 
terval) . 

We discussed in App.lXlthat a truncated, finite Cheby- 
shev expansion of Ab(u>) may fail to be positive although 
ylg(w) is. Positivity of the finite Chebyshev expansion 
is ensured by using attenuation factors of a positive ker- 
nel like in Eq. (|A5|) . To obtain a positive definite scalar 
product (|C3[) for a finite Chebyshev expansion we must 
therefore use modified moments /i^ g„ instead of unmod- 
ified moments /j,^. For the Lanczos algorithm, used in 
Sec. I Villi which aims at constructing an orthonormal 
basis in Krylov-subspaces, positive definiteness of (|C3|) 
is crucial. For the calculation of spectral or dynamical 
properties it is not that essential. But we learned in 
App-E]that the modification of moments by attenuation 
factors never ruins a result. Hence our recommendation 
is to always modify the moments [i„ put into the calcu- 
lation with attenuation factors. 



APPENDIX D: EIGENSTATES OF H% 



term clm-i\M) is missing, and the scalar product is not 
symmetric in | a) , | b) . The mathematical reason that Hg 
fails to be hermitian, is that Hg is obtained from Hb 
via a non-orthogonal projection. Working with a non- 
hermitian operator could in principle spoil the calculation 
for M < N. We found in Sec. [V]that this does not hap- 
pen. By an explicit calculation of the spectral function 
encoded by Hg , we now show that the non-hermiticity 
of Hg has no (negative) consequences. 

As Fig. [5] shows, for N > M CSM resolves M poles in 
the spectral function encoded by Hg . The position of 
the poles is determined by Xj , and their weight Wj can be 
calculated using the eigenstates \(f)j). We first expand the 
th Chebyshev vector |0) in the \(j>j). With the identity 



1 M \ - 8; 

y ^, T m (xj)T n (xj) = < ^ 



to, n 
to = n = ' 



(D3) 

which is a kind of discrete orthogonality relation for 
Chebyshev polynomials, we find that 



M 

I°> = m5>>- 



(D4) 



j=i 



The coefficients in this linear combination are To (xj ) = 1 . 

Now S(l> — Hg)\<f)j) = 6(u> — Xj)\<j)j), and we obtain 

M 
3=1 



According to App. [B] the characteristic polyno- 
mial of Hg 1 (Eq. J26])) is det[x - [H^) rnn ] = 
2-( M - 1 )T M (x). The M th Chebyshev polynomial 
Tm{x) = cos(Marccosa;) has M distinct real roots 



1.. 



(Dl) 



We concluded in Sec. [V] that Hg is diagonalizable with 
real eigenvalues Xj . The eigenstates of Hg can be given 
explicitly, using the recurrence ([T]). The state 



l^> = 10) 



M-l 

! E 

m— 1 



T m (xj)\m) 



is the eigenstate of H B to eigenvalue Xj, i.e 
Xj\4>j). If we use Eq. 



3 ' 

to calculate -ffjf | 



^1 



we find 



that the term proportional to Tm{Xj) drops out, because 
Xj is a root of Tm{x). This explains why the roots of 
Tm{x) occur as eigenvalues of Hg . 

Scalar products of the states \(j>j) have to be calculated 
using Eq. (|C3p . It turns out, that the \<pj) are not or- 
thogonal to each other. Indeed, while Hb is hermitian, 
the truncated Hg is not hermitian. If one repeats the 
calculations leading to Eq. (|C5|) . now with upper sum- 
mation boundary M — 1 instead of oo, one finds that a 



with the weight of each pole given by 



M 



<0|^> = 



1 

M 



M—l 



2 T m(Xj)fJ>r, 
m—l 



(D6) 



The first line in this equation follows from the defini- 
tion (|D2p , the second line from comparison with Eq. (0| . 
Here Ab(lu) is assumed as an expansion with M Cheby- 
shev moments, in accordance with the truncation of Hg . 
The weight Wj of the pole at Xj is therefore the value of 
Ab{u) at the position of the pole, weighted by the in- 
(D2) verse 

{l-x 2 ) 1/2 of the weighting function for Chebyshev 
expansions (and with 1/M for M poles). Using the at- 
tenuation factors from Eq. (|A5|) . Ag(w) is positive, which 
implies that the weight Wj of each pole is positive. The 
total weight is 



M M 



(D7) 



according to the sum rule /Iq 3 = f_~ Ab{x)cIx = 1. This 
result also follows with the explicit expression for Wj in 
terms of the n% (first line of Eq. JD6|), and Eq. (|D3]l . If 



we express Wj by Ab(u>) (second line of Eq. (|D6|), we find 
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a nice relation to quadrature formulas from numerical 
integration. In the expression for the total weight 

M M 

E^-^^EV 1 "^ 2 ^^ ' (D8) 

the right hand side is the formula for Gauss-Chebyshev 
integration 2 ^ of a function, with abscissas at Xj. Now 
Ag(u), expanded with M Chebyshev moments, is a 
polynomial of order M times the weighting function 
(1 — a; 2 ) -1 / 2 . The Gauss-Chebyshev integration formula 
is exact in this case, and we find once again J2j=i w j 
J Ab{x)cLx = 1. 

We can finally understand why if|f is non-hermitian 
using the scalar product Eq. (|C9|) . A state \<j>j) corre- 
sponds to the function 

M-l 

) (D9) 

m— 1 

which is the M th Chebyshev approximation to 8(x — Xj). 
With a finite number of moments, (l)j( x ) has finite width. 
Using Eq. (|C9|) for the scalar product, 

{(j )j \4> k )=J (f> j {x)<f>k(x)A B (x)dx , (D10) 



we see that the integral is non-zero even for j ^ k, since 
the two functions (f>j(x), 4>k{x) have finite overlap. The 
deviation of Hg from hermiticity is therefore an effect 
of finite resolution of a finite Chebyshev expansion. The 
larger M, the smaller is the overlap and hence the scalar 
product. We conclude that in the limit M — > oo, the 
states \4>j) are mutually orthogonal, in accordance with 
hermiticity of Hb proved in App. [C] 



We have obtained two important results in this Ap- 
pendix. First, despite its non- hermiticity, Hg encodes 
a positive, normalized spectral function. Using the trun- 
cated Hg, as in Sec. |V] for M < N, does therefore not 
lead to erroneous results. Second, we related the CS con- 
struction for Hb to an explicit representation Eq. ([7]) . In 
particular, Eq. (|D6j) shows how to obtain from given mo- 
ments tt^j a discretization (Eq. (fT2)l ) of Ab{uj) (the op- 
posite way is given by Eq. (O). The point is that the CS 
construction indicates how to discretize a given A^(w) to 
prescribed resolution. Within CSM we can work without 
a discretization anyway. The reader should note that the 
CS construction is no longer equivalent to the discretiza- 
tion of Ab{oj) in the extensions of CSM mentioned in the 
Summary (Sec. |X|). 
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